{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0a108f97",
   "metadata": {},
   "source": [
    "\n",
    "# Functional Bayesian Networks\n",
    "\n",
    "Functional Bayesian Networks (FBNs) are Bayesian networks where each CPD is a Python function that returns a Pyro distribution. This lets you model arbitrary discrete, continuous, or mixed relationships—while keeping standard graph semantics for sampling, interventions, and learning.\n",
    "\n",
    "Similar to other Bayesian Network classes in pgmpy, there are two main components of the model: 1. The graphical structure of the model, 2. The parameterization of the model (defined using `FunctionalCPD`).\n",
    "\n",
    "This tutorial introduces **Functional Bayesian Networks (FBNs)** and the accompanying **Functional CPDs** in `pgmpy`. You'll learn how to:\n",
    "- Define FunctionalCPDs as *functions that return Pyro distributions*.\n",
    "- Build FunctionalBayesianNetworks on mixed data.\n",
    "- Simulate data from the model.\n",
    "- Use vectorized CPDs for performance.\n",
    "- Fit simple parametric FBNs using **SVI** and **MCMC** via Pyro."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "164702d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import torch\n",
    "import pyro\n",
    "import pyro.distributions as dist\n",
    "\n",
    "from pgmpy.models import FunctionalBayesianNetwork\n",
    "from pgmpy.factors.hybrid import FunctionalCPD\n",
    "\n",
    "# Reproducibility\n",
    "pyro.set_rng_seed(123)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0455f441",
   "metadata": {},
   "source": [
    "\n",
    "## Functional CPDs: the core idea\n",
    "\n",
    "A **Functional CPD** has three arguments:\n",
    "1. `variable`: The variable for the FunctionalCPD is being defined. For a fully parameterized model, each node in the model needs a `FunctionalCPD` associated with it.\n",
    "2. `fn`: A Python callable that takes a dictionary of parents' values as input and returns a **Pyro distribution**.\n",
    "3. `parents`: The parents of `variable` in the model.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "215bc744",
   "metadata": {},
   "source": [
    "## Example 1: A simple Gaussian chain: `x1 → x2 → x3`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2b606340",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "      <th>x3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.111467</td>\n",
       "      <td>1.015461</td>\n",
       "      <td>1.525792</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.120363</td>\n",
       "      <td>0.610113</td>\n",
       "      <td>0.519270</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.369635</td>\n",
       "      <td>0.326770</td>\n",
       "      <td>0.861743</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.240418</td>\n",
       "      <td>0.969617</td>\n",
       "      <td>1.934877</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-1.196924</td>\n",
       "      <td>-0.011801</td>\n",
       "      <td>0.641020</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         x1        x2        x3\n",
       "0 -0.111467  1.015461  1.525792\n",
       "1  0.120363  0.610113  0.519270\n",
       "2 -0.369635  0.326770  0.861743\n",
       "3 -0.240418  0.969617  1.934877\n",
       "4 -1.196924 -0.011801  0.641020"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "# Build the structure\n",
    "gauss_chain = FunctionalBayesianNetwork([(\"x1\", \"x2\"), (\"x2\", \"x3\")])\n",
    "\n",
    "# Define CPDs\n",
    "cpd_x1 = FunctionalCPD(\"x1\", fn=lambda _: dist.Normal(0.0, 1.0))  # prior\n",
    "cpd_x2 = FunctionalCPD(\"x2\", fn=lambda parents: dist.Normal(1.0 + 0.8 * parents[\"x1\"], 0.5), parents=[\"x1\"])\n",
    "cpd_x3 = FunctionalCPD(\"x3\", fn=lambda parents: dist.Normal(0.3 + 1.0 * parents[\"x2\"], 1.0), parents=[\"x2\"])\n",
    "\n",
    "gauss_chain.add_cpds(cpd_x1, cpd_x2, cpd_x3)\n",
    "gauss_chain.check_model()\n",
    "\n",
    "# Draw a few samples\n",
    "samples_gc = gauss_chain.simulate(n_samples=5, seed=123)\n",
    "samples_gc\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a727558d",
   "metadata": {},
   "source": [
    "### Example 2: Complex model with mixture data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bb1b52e5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "      <th>w</th>\n",
       "      <th>y</th>\n",
       "      <th>z</th>\n",
       "      <th>c</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.111467</td>\n",
       "      <td>0.888683</td>\n",
       "      <td>-0.363940</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-0.126706</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.120363</td>\n",
       "      <td>0.369773</td>\n",
       "      <td>-0.469728</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.203071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.369635</td>\n",
       "      <td>0.752397</td>\n",
       "      <td>-0.719913</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.660035</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.240418</td>\n",
       "      <td>0.030989</td>\n",
       "      <td>-0.391040</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.576696</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-1.196924</td>\n",
       "      <td>0.781968</td>\n",
       "      <td>-1.086602</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.384325</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.209269</td>\n",
       "      <td>1.298313</td>\n",
       "      <td>0.468003</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.734681</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>-0.972355</td>\n",
       "      <td>0.923385</td>\n",
       "      <td>-1.151972</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-0.902663</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>-0.755045</td>\n",
       "      <td>1.667385</td>\n",
       "      <td>-1.473914</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-0.654092</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         x1        x2         w    y    z         c\n",
       "0 -0.111467  0.888683 -0.363940  0.0  1.0 -0.126706\n",
       "1  0.120363  0.369773 -0.469728  0.0  0.0  0.203071\n",
       "2 -0.369635  0.752397 -0.719913  0.0  0.0  0.660035\n",
       "3 -0.240418  0.030989 -0.391040  0.0  1.0  0.576696\n",
       "4 -1.196924  0.781968 -1.086602  0.0  0.0  0.384325\n",
       "5  0.209269  1.298313  0.468003  1.0  0.0  1.734681\n",
       "6 -0.972355  0.923385 -1.151972  0.0  1.0 -0.902663\n",
       "7 -0.755045  1.667385 -1.473914  1.0  1.0 -0.654092"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "complex_bn = FunctionalBayesianNetwork(\n",
    "    [(\"x1\", \"w\"), (\"x2\", \"w\"), (\"x1\", \"y\"),\n",
    "     (\"x2\", \"y\"), (\"w\", \"y\"), (\"y\", \"z\"), \n",
    "     (\"w\", \"z\"), (\"y\", \"c\"), (\"w\", \"c\")]\n",
    ")\n",
    "\n",
    "# Roots\n",
    "cpd_x1 = FunctionalCPD(\"x1\", fn=lambda _: dist.Normal(0.0, 1.0))\n",
    "cpd_x2 = FunctionalCPD(\"x2\", fn=lambda _: dist.Normal(0.5, 1.2))\n",
    "\n",
    "# Continuous mediator: w = 0.7*x1 - 0.3*x2 + ε\n",
    "cpd_w = FunctionalCPD(\n",
    "    \"w\",\n",
    "    fn=lambda parents: dist.Normal(0.7 * parents[\"x1\"] - 0.3 * parents[\"x2\"], 0.5),\n",
    "    parents=[\"x1\", \"x2\"]\n",
    ")\n",
    "\n",
    "# Bernoulli target with logistic link: y ~ Bernoulli(sigmoid(-0.7 + 1.5*x1 + 0.8*x2 + 1.2*w))\n",
    "cpd_y = FunctionalCPD(\n",
    "    \"y\",\n",
    "    fn=lambda parents: dist.Bernoulli(logits=(-0.7 + 1.5 * parents[\"x1\"] + 0.8 * parents[\"x2\"] + 1.2 * parents[\"w\"])),\n",
    "    parents=[\"x1\", \"x2\", \"w\"]\n",
    ")\n",
    "\n",
    "# Downstream Bernoulli influenced by y and w\n",
    "cpd_z = FunctionalCPD(\n",
    "    \"z\",\n",
    "    fn=lambda parents: dist.Bernoulli(logits=(-1.2 + 0.8 * parents[\"y\"] + 0.2 * parents[\"w\"])),\n",
    "    parents=[\"y\", \"w\"]\n",
    ")\n",
    "\n",
    "# Continuous outcome depending on y and w: c = 0.2 + 0.5*y + 0.3*w + ε\n",
    "cpd_c = FunctionalCPD(\n",
    "    \"c\",\n",
    "    fn=lambda parents: dist.Normal(0.2 + 0.5 * parents[\"y\"] + 0.3 * parents[\"w\"], 0.7),\n",
    "    parents=[\"y\", \"w\"]\n",
    ")\n",
    "\n",
    "complex_bn.add_cpds(cpd_x1, cpd_x2, cpd_w, cpd_y, cpd_z, cpd_c)\n",
    "complex_bn.check_model()\n",
    "\n",
    "# Simulate data from it\n",
    "complex_bn.simulate(n_samples=8, seed=123)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8329345",
   "metadata": {},
   "source": [
    "\n",
    "## Vectorized CPDs for speed\n",
    "\n",
    "Set `vectorized=True` and have your `fn(parent_df)` return a **batched** Pyro distribution whose batch size equals the number of rows in the provided parent DataFrame. This makes sampling much faster for large `n_samples`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c65e71cb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.111467</td>\n",
       "      <td>0.796616</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.120363</td>\n",
       "      <td>0.786959</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.369635</td>\n",
       "      <td>-0.077290</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.240418</td>\n",
       "      <td>0.097413</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-1.196924</td>\n",
       "      <td>-0.559637</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         x1        x2\n",
       "0 -0.111467  0.796616\n",
       "1  0.120363  0.786959\n",
       "2 -0.369635 -0.077290\n",
       "3 -0.240418  0.097413\n",
       "4 -1.196924 -0.559637"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "from pgmpy import config\n",
    "\n",
    "vec_bn = FunctionalBayesianNetwork([(\"x1\", \"x2\")])\n",
    "\n",
    "cpd_x1 = FunctionalCPD(\"x1\", fn=lambda _: dist.Normal(0.0, 1.0))\n",
    "\n",
    "def x2_fn_vec(P):\n",
    "    x1 = torch.tensor(P[\"x1\"].values, dtype=config.get_dtype(), device=config.get_device())\n",
    "    mu = 0.5 + 0.9 * x1\n",
    "    sigma = torch.full_like(mu, 0.3)\n",
    "    return dist.Normal(mu, sigma)\n",
    "\n",
    "cpd_x2 = FunctionalCPD(\"x2\", fn=x2_fn_vec, parents=[\"x1\"], vectorized=True)\n",
    "\n",
    "vec_bn.add_cpds(cpd_x1, cpd_x2)\n",
    "vec_bn.check_model()\n",
    "\n",
    "# Large draw to highlight performance of vectorized CPDs\n",
    "vec_samples = vec_bn.simulate(n_samples=5000, seed=123)\n",
    "vec_samples.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be69122b",
   "metadata": {},
   "source": [
    "\n",
    "## Parameter learning with SVI\n",
    "\n",
    "When CPDs contain **Pyro parameters** (`pyro.param(...)`), you can fit them to data using `model.fit(..., method=\"SVI\")`.\n",
    "Below, we synthesize data from a simple linear-Gaussian model and then recover the parameters.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c92c26bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:pgmpy:Step 0 | Loss: 6574.2528\n",
      "INFO:pgmpy:Step 50 | Loss: 4983.4878\n",
      "INFO:pgmpy:Step 100 | Loss: 4229.9033\n",
      "INFO:pgmpy:Step 150 | Loss: 3987.9143\n",
      "INFO:pgmpy:Step 200 | Loss: 3981.4799\n",
      "INFO:pgmpy:Step 250 | Loss: 3981.4706\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'x1_mu': 0.7751224637031555,\n",
       " 'x1_sigma': 0.6131476759910583,\n",
       " 'x2_inter': 1.1970738172531128,\n",
       " 'x2_sigma': 0.6990763545036316}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Generate synthetic data\n",
    "true_mu, true_sigma = 0.8, 0.6\n",
    "N = 2000\n",
    "x1 = torch.normal(mean=true_mu, std=true_sigma, size=(N,))\n",
    "# FIX: vectorized draw for x2 (no size when mean is a tensor)\n",
    "x2 = torch.normal(mean=1.2 + x1, std=0.7)         # or: x2 = 1.2 + x1 + 0.7 * torch.randn_like(x1)\n",
    "\n",
    "data = pd.DataFrame({\"x1\": x1.numpy(), \"x2\": x2.numpy()})\n",
    "\n",
    "from torch.distributions import constraints\n",
    "import pyro, pyro.distributions as dist\n",
    "pyro.clear_param_store()  # helpful if you re-run the cell\n",
    "\n",
    "svi_bn = FunctionalBayesianNetwork([(\"x1\", \"x2\")])\n",
    "\n",
    "def x1_fn(_):\n",
    "    mu = pyro.param(\"x1_mu\", torch.tensor(0.0))\n",
    "    sigma = pyro.param(\"x1_sigma\", torch.tensor(1.0), constraint=constraints.positive)\n",
    "    return dist.Normal(mu, sigma)\n",
    "\n",
    "def x2_fn(p):\n",
    "    inter = pyro.param(\"x2_inter\", torch.tensor(0.0))\n",
    "    sigma = pyro.param(\"x2_sigma\", torch.tensor(1.0), constraint=constraints.positive)\n",
    "    return dist.Normal(inter + p[\"x1\"], sigma)\n",
    "\n",
    "svi_bn.add_cpds(\n",
    "    FunctionalCPD(\"x1\", fn=x1_fn),\n",
    "    FunctionalCPD(\"x2\", fn=x2_fn, parents=[\"x1\"]),\n",
    ")\n",
    "svi_bn.check_model()\n",
    "\n",
    "# Fit with SVI\n",
    "params_svi = svi_bn.fit(data, method=\"SVI\", num_steps=300)\n",
    "{k: v.item() if torch.is_tensor(v) and v.ndim == 0 else v for k, v in params_svi.items()}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c20f73f",
   "metadata": {},
   "source": [
    "\n",
    "## Bayesian inference with MCMC\n",
    "\n",
    "If you prefer a fully Bayesian approach, provide **priors** and switch `method=\"MCMC\"`.  \n",
    "Here, we place priors on the parameters, then sample from the posterior.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "41483dd0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sample: 100%|█████████████████████████████████████████| 400/400 [00:03, 109.04it/s, step size=6.89e-01, acc. prob=0.881]                                                                                                                                                                            \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'x1_mu': (0.777364034922299, 0.0144704934171712),\n",
       " 'x1_sigma': (0.6129607990640703, 0.010122547469210461),\n",
       " 'x2_inter': (1.1970455001675568, 0.014884526650204519),\n",
       " 'x2_sigma': (0.7008062278971647, 0.010433891834526142)}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mcmc_bn = FunctionalBayesianNetwork([(\"x1\", \"x2\")])\n",
    "\n",
    "# Priors with matched dtype/device (avoid dtype mismatches)\n",
    "dtype = config.get_dtype()\n",
    "device = config.get_device()\n",
    "def prior_fn():\n",
    "    t = lambda v: torch.tensor(v, dtype=dtype, device=device)\n",
    "    return {\n",
    "        \"x1_mu\":   dist.Normal(t(0.0), t(5.0)),\n",
    "        \"x1_sigma\": dist.HalfNormal(t(2.0)),\n",
    "        \"x2_inter\": dist.Normal(t(0.0), t(5.0)),\n",
    "        \"x2_sigma\": dist.HalfNormal(t(2.0)),\n",
    "    }\n",
    "\n",
    "# CPDs consume *sampled* prior values (from the model) + parents\n",
    "def x1_fn_prior(priors, _):\n",
    "    return dist.Normal(priors[\"x1_mu\"], priors[\"x1_sigma\"])\n",
    "\n",
    "def x2_fn_prior(priors, p):\n",
    "    return dist.Normal(priors[\"x2_inter\"] + p[\"x1\"], priors[\"x2_sigma\"])\n",
    "\n",
    "mcmc_bn.add_cpds(\n",
    "    FunctionalCPD(\"x1\", fn=x1_fn_prior),\n",
    "    FunctionalCPD(\"x2\", fn=x2_fn_prior, parents=[\"x1\"]),\n",
    ")\n",
    "\n",
    "pyro.clear_param_store()\n",
    "\n",
    "post = mcmc_bn.fit(\n",
    "    data,\n",
    "    method=\"MCMC\",\n",
    "    prior_fn=prior_fn,\n",
    "    num_steps=200,\n",
    "    nuts_kwargs={\"target_accept_prob\": 0.8},\n",
    "    mcmc_kwargs={\"num_chains\": 1, \"warmup_steps\": 200},\n",
    ")\n",
    "\n",
    "# Peek at posterior summaries\n",
    "{k: (v.mean().item(), v.std().item()) for k, v in post.items() if torch.is_tensor(v)}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9572436f",
   "metadata": {},
   "source": [
    "\n",
    "## Interventions and conditioning (preview)\n",
    "\n",
    "Conceptually, an intervention `do(X = value)` **severs incoming edges into `X`** and replaces its CPD with a constant or a new distribution.  \n",
    "A simple approach for simulations is to *temporarily swap the CPD of `X`* with a `FunctionalCPD` that ignores parents and returns `dist.Delta(value)` (or any desired distribution).\n",
    "\n",
    "Conditioning on continuous **point evidence** should use **likelihood weighting** or proper inference (SVI/MCMC), not plain rejection sampling. The high-level recipe for likelihood weighting is:\n",
    "1. Sample all **non-evidence** nodes in topological order (respecting any `do(...)` replacements).\n",
    "2. Clamp evidence nodes to their observed values.\n",
    "3. Weight each draw by the product of evidence likelihoods under their parents’ sampled values.\n",
    "4. Normalize weights and either keep **weighted samples** or **resample** for an unweighted posterior sample.\n",
    "\n",
    "Future versions may expose `simulate(do=..., evidence=...)` directly in the API.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d04cfb3f",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "### Key takeaways\n",
    "- Functional CPDs let you specify **any** distribution you can write in Pyro.\n",
    "- Mixed types (discrete/continuous) are straightforward.\n",
    "- Use **vectorized CPDs** for performance on large simulations.\n",
    "- For learning: quick **SVI** with `pyro.param(...)` or fully Bayesian **MCMC** with `prior_fn`.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
